Nature of the quantum phase transitions in the two-dimensional hardcore boson model 
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We use two Quantum Monte Carlo algorithms to map out the phase diagram of the two- 
dimensional hardcore boson Hubbard model with near (Vi) and next near (V2) neighbor repulsion. 
At half filling we find three phases: Superfluid (SF), checkerboard solid and striped solid depend- 
ing on the relative values of Vi, V2 and the kinetic energy. Doping away from half filling, the 
checkerboard solid undergoes phase separation: The superfluid and solid phases co-exist but not as 
a single thermodynamic phase. As a function of doping, the transition from the checkerboard solid 
is therefore first order. In contrast, doping the striped solid away from half filling instead produces 
a striped supersolid phase: Co-existence of density order with superfluidity as a single phase. One 
surprising result is that the entire line of transitions between the SF and checkerboard solid phases 
at half filling appears to exhibit dynamical 0(3) symmetry restoration. The transitions appear 
to be in the same universality class as the special Heisenberg point even though this symmetry is 
explicitly broken by the V2 interaction. 

PACS numbers: 75.10Nr, 05.30 Jp, 67.40.Yv, 74.60.Ge 



I. INTRODUCTION 



The boson Hubbard HamiltonianEJ has been studied 
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as a model of the superconductor-insul 
in materials with preformed Cooper pa^s, 
in disordered and restricted geometries 
transitions in quantum spin systems in <p 
netic fieldSjEil and of supersolid behaviorE 
the fermion Hubbard Hamiltonian, the boson model ex- 
plores the role of correlations in inducing ordered phases 
of many quantum mechanical particles, and the nature 
of the quantum phase transitions between these phases. 
However, unlike the fermion case where it is very dif- 
ficult to reach low temperatures away from points of 
special particle-hole symmetry, Quantum Monte Carlo 
(QMC) simulations of the doped boson system have no 
"sign problem" and hence can successfully be performed. 
Many fascinating and unexpected features pjise, for ex- 
ample re-entojijL-Mott insulating b ehayi pio, universal 
conductivityJilclEa and supersolidityOtL^I Indeed, until 
algorithms are developed to deal with the sign problem in 
fermion QMC the boson Hubbard Hamiltonian offers us 
the best opportunity to explore systematically the details 
of the competition between phases with diagonal and off- 
diagonal long range order with QMC simulations. 

In this paper we extend some of the previous work 
which- established the basic phase diagram of the 
mode]o't£rE£l in order to characterize the detailed criti- 
cal properties of the transitions between the phases. We 



will use two recently developed algorithms. The firsttSEHl 
is based on a duality transformation which enables an ex- 
act mapping of the boson Hubbard Hamiltajjiaji onto a 
model of conserved currents, and the secondE3E3 is based 
on a stochastic series expansion of the imaginary time 
evolution operator. 

The paper is organized as follows: We will first review 
the definition of the boson Hubbard model and some of 
its basic qualitative properties. We will then provide a 
brief review of the two numerical algorithms we employ. 
Our presentation of the results begins with a discussion of 
the half-filled system which focusses on a scaling analysis 
of the transition from a strong coupling solid to a weak 
coupling superfluid phase. Next, away from half-filling, 
we discuss the nature of the coexistence of the solid and 
the superfluid. Finally, we present some concluding re- 
marks and open questions. 



II. THE BOSON HUBBARD MODEL 

The hardcore boson Hubbard Hamiltonian is, 

(iJ> 

+ Vi njUj + Vi mnk- (1) 

(U> «i,k» 
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cti (oT ) are destruction (creation) operators of hard-core 
bosons on site i of a 2-d L x L square lattice, and m is 
the density at site i. The hopping parameter is chosen to 
be t — 1 to fix the energy scale. V\ (V2) is the near (next 
near) neighbor interaction. 

For weak couplings, the ground state of the Hamil- 
tonian is a supcrfluid. Increasing the near neighbor in- 
teraction strength, V\, at half filling drives a transition 
to a checkerboard solid phase where the sites are alter- 
nately occupied and empty. This phase is characterized 
by a vanishing superfluid density, p s , long range density- 
density correlations, and a gap in the energy spectrum 
reflected, for example, as a vanishing compressibility, k, 
and corresponding plateau in a plot of density p versus 
chemical potential /i. Increasing the next near neighbor 
interaction strength V2 at half filling can likewise drive 
a transition to a striped solid, where horizontal (or ver- 
tical) lines of sites are alternately occupied. This phase 
also has p s — 0, k — and long range density-density 
correlations. 

At V2 = 0, and after an appropriate sublattice spin 
rotation, the hard-core boson model is equivalent to the 
spin-i antiferromagnetic XXZ model, with half-filling 
corresponding to the zero magnetization sector. The hop- 
ping t maps onto J x /2 while the interaction strength V\ 
maps onto J z and the chemical potential p is related 
to the magnetic field as h = p — zVi/2, where z is the 
number of nearest neighbors of a lattice site. In this lan- 
guage, superfluid order corresponds to magnetic order in 
the XY plane, while density order corresponds to mag- 
netic order in the Z direction. The boson superfluid-cdw 
insulator phase transition at V\ = 2t corresponds to the 
XY-Ising change in universality class at J x — J z . At this, 
the Heisenberg point, the Hamiltonian has an 0(3) sym- 
metry and as a consequence the critical temperature is 
driven to T c = 0. This symmetry is explicitly broken for 
other values of V\ or for nonzero V%. The boson-Hubbard 
model with nonzero V2 also has a spin analog, namely to 
a Hamiltonian with next-near neighbor exchange. 

The behavior of the boson-Hubbard model away from 
half-filling is considerably more complex. While the com- 
pressibility surely becomes nonzero, so that the state is 
not a Mott insulator, it is possible that, despite doping, 
the charge correlations remain long ranged. If the doped 
holes are mobile and interspersed with the density or- 
dered bosons, one has simultaneous superfluid order. On 
the other hand, the doped holes might phase separate 
leaving distinct regions of the lattice with superfluid and 
charge ordering. We will carefully explore these alternate 
possibilities and describe the nature of the transitions be- 
tween them. 



III. THE MONTE CARLO ALGORITHMS 

Until relatively recently, algorithms to simulate inter- 
acting quantum bosons on a lattice suffered significant 



weaknesses. Most importantly, they had problems with 
extremely long autocorrelation times, which, as in clas- 
sical Monte Carlo, were caused by the inability of lo- 
cal changes to move configurations effectively through 
phase space. A related difficulty was that local moves 
also resulted in global conservation laws which limited 
the accessible regions of phase space and hence the mea- 
surements that could be performed. In the last several 
years, "cluster" and "loop" algorithms have very success- 
fully addressed some of these problems. Efl Unfortunately, 
these algorithms do not work equally well in all regions of 
parameter space, notably away from half filling. In addi- 
tion, they are not simple to implement when the Hamil- 
tonian is made more complex, e.g. with the inclusion of 
disorder or longer range interactions.E3 In this section we 
review two approaches which have short autocorrelation 
times but also are easily implemented for these more gen- 
eral models. 

Dual Algorithm: To perform efficient boson simula- 
tions we use a newly developed QMC algorithm based 
on the exact duality transformation of the Bosonic Hub- 
bard modelOtJ This approach begins by expressing the 
partition function as a path integral over coherent states. 
To implement the duality transformation exactly we fol- 
lowed the method of Ref. |2f] The result is that hard- 
core bosons are represented by conserved "currents" that 
propagate in the positive imaginary time direction and 
which can make jumps in any spatial direction. The path 
integral is transformed into a sum over all deformations 
of the currents, very much like the world line algorithm. 
This formulation is similar to that of Refs. |^,|| in that it 
relies on the duality transformation, but different in that 
this transformation is exact in our case. Consequently, 
with our algorithm we simulate the true hard-core Hub- 
bard model, whereas Refs. |^|| study the very high density 
(many bosons/site) limit. 

The simulation is done in the standard way. A defor- 
mation is suggested and it is accepted or rejected in a 
way satisfying detailed balance. A feature of this algo- 
rithm which is not easily implemented in the world line 
algorithm is that the imaginary time step is not subdi- 
vided by a checkerboard break-up, and so a particle can 
hop several lattice sites at a time. This has the great 
advantage that it enables us to measure the correlation 
function of the superfluid order parameter, (ajaj), for 
large |i — j|, which is very difficult to measure efficiently 
in older approaches. 

Stochastic Series Expansion: The stochastic series-epf-i 
pansion (SSE) algorithm with operator loop updateC3'E2l 
was used for additional calculations on larger systems in 
the grand canonical ensemble. This algorithm does not 
suffer from time discretization errors and is one of the 
most effective algorithm for quantum systems. It has 
been shown to be as effective as the loop algorithrrO 
for models where the loop algorithm can be applied. In 
bosonic models it is better than the loop algorithm, since 
it does not suffer from the exponential slowing down of 
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the loop algorithm when the chemical potential is tuned 
away from half filling.E3. 

Measurements: The energy E is obtained as the expec- 
tation value of the Hamiltonian and is used to determine 
the chemical potential via the relation p = dE/dN, from 
simulations in the canonical ensemble. Hysteretic behav- 
ior in E also provides supporting evidence in character- 
izing the order of phase transitions. 

It is straightforward to measure the density-density 
correlations and their Fourier transform, the structure 
factor: 

c(l) = 

S(q) = ^e^ 1 c(l)- (2) 
1 

These quantities characterize the diagonal long range or- 
der. Divergence of S(ir, tt) in the thermodynamic limit 
indicates checkerboard order, and in S(jr,0) or 5(0, 7r) 
striped order. On finite lattices, the structure factor at 
the appropriate momentum diverges as the system size 
in the ordered phase, so that a scaling analysis of simula- 
tions on different lattice sizes can demonstrate long range 
order. From the density itself comes the compressibility 
which characterizes Mott insulating behavior. 

It is also crucial to obtain the winding number, since 
it will be used to determine the superfluid density p s oc 
(W 2 ), a relationship first emphasized in the context 
of Quantum Monte Carlo simulations by Ceperley and 
Pollock.EEl While global conservation laws on winding 
and particle numbers preclude the straightforward evalu- 
ation of the superfluid density in the dual algorithm and 
other traditional world-line approaches, we used meth- 
ods which circumvent this difficultyJl3't3 Specifically, we 
calculate the "pseudo" current-current correlation func- 
tion at different imaginary times. The "pseudocurrent" 
is defined as the net number of bosons which jumped 
in a given direction in one imaginary time step. It can 
be shown easily that the Fourier transform of this cor- 
relation function approac hcSpXVF 2 ) as the (Matsubara) 
frequency approaches zero.E£Ej 

Meanwhile, the global loop updates of the SSE algo- 
rithm allow modifications of the particle number and of 
the winding numbers. The superfluid density p s can thus 
be measured directly from the winding number fluctua- 
tions in these approaches. Two of us have recently devel- 
oped a method to measure Green's functions like (<2jOj) 
in the SSE algorithm. We refer to Ref. ^3| for details of 
this method. Values for the superfluid density obtained 
via the dual and SSE algorithms, although obtained very 
differently, are entirely consistent. 

IV. PHASE DIAGRAM AND TRANSITIONS 
HALF FILLING 

As we have already described, at half filling, one can 
easily recognize the existence of at least three phases. For 
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FIG. 1. The structure factor, S(tt,tt) (circles), S(n,0) 
(squares) and the superfluid density, p a (triangles) at half 
filling for Vi = 6 as functions of V 2 . L = 12, /? = 12 

weak V\ and V2 it is clear that a superfluid phase exists 
and, as we will see below, this extends to very strong re- 
pulsions when the competing interactions nearly balance. 
If Vi dominates, the energy cost of having near neighbors 
becomes too high and the bosons organize themselves 
into a checkerboard solid. On the other hand, when V2 
dominates, it is less costly to have near neighbors com- 
pared to next near neighbors and the bosons organize 
themselves in a striped solid. What is not obvious is 
whether there are other phases, for example supersolids 
separating the, solid phases from the superfluid phase. In 
previous workE3E3E3 we demonstrated that at half fill- 
ing there are no supersolid phases in this model. This 
is confirmed in the present work where we determine the 
phase diagram more accurately than before and study in 
detail the nature of the transitions. The case V2 = and 
V\ = 2 serves as a good test for our simulations, since it 
corresponds to the well-understood Heisenberg point of 
the XXZ model. 

The checkerboard and striped solids have structure fac- 
tors which diverge as L 2 , with momentum ordering vec- 
tors (71-, 71") and (tt, 0) (or (0,7r)) respectively. Figure 1 
shows these two quantities for fixed V\ as V2 is varied. We 
see that as V2 is increased, 5(7r,7r) falls to zero, p s takes 
a finite value while 5(7r, 0) remains zero. This indicates a 
(tt, 7r)-solid to superfluid phase transition. Increasing V2 
further, p s vanishes abruptly while S(tt, 0) takes a finite 
value indicating a superfluid to striped (tt, 0)-solid. 

Putting together the transition points obtained from 
several such slices, we arrive at the ground state phase 
diagram in the (Vi , V2) plane at half filling. This is shown 
in Fig. 2. Two remarks are in order. It is interesting to 
note in Fig. 2 that even for very large values of V\ and 
V2 , there is no direct transition between the checkerboard 
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FIG. 2. The ground state phase diagram in the (Vi,V2) 
plane. The dashed line indicates first order transitions, the 
solid line exhibits what appears to be dynamically restored 
0(3) symmetry except at (Vi = 2, V2 = 0) where the symme- 
try is explicit in the hamiltonian (see text below). 



and striped phases: The superfluid phase seems always to 
intervene. Presumably, the two solid phases eventually 
meet as both Vi and V% gQ-to infinity. This is in contrast 
with the mean field resultEj that the two phases meet at 
(Vi, V2) = (4, 2). The second remark is that V2 does not 
need to be larger than Vi in order for striped order to win 
over checkerboard order. For example, (Vi, V2) = (6, 4) is 
in the striped phase. This is important because it makes 
the striped phase more likely to appear physically since 
one might expect near neighbor repulsion to be stronger 
than next near neighbor. 

Fig. 1 also shows that (ir, 7r)-SF transition appears 
smooth suggesting a continuous phase transition. This 
will be examined in detail below. On the other hand, 
Fig. 1 also shows that the SF-(7r, 0) transition is very 
sudden, suggesting a first order transtion. To verify the 
first order nature of the SF-(-7r, 0) transition we look for 
hysteresis in S(ir, 0) and (E) as V2 is increased and then 
decreased, always starting a new simulation from the 
last configuration of the previous run. The results are 
shown in Fig. 3. Hysteresis is clearly seen, supporting 
the evidence for a first order transition. Similar param- 
eter sweeps through the (it, 7r)-SF transition exhibit no 
hysteresis. 

To determine the order of the transition from checker- 
board solid to superfluid we do finite size scaling. Since 
the Heisenberg point (Vi = 2, V 2 = 0) is very special (ex- 
plicit 0(3) in the Hamiltonian), we first did the analysis 
away from it by fixing Vi = 3 and finding the transition 
as V2 is changed (see Fig. 1). Since from Fig. 1, the SF 
to (7r, 7r)-solid transition appears continuous, we will first 
proceed by assuming a second order transition and car- 
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FIG. 3. Hysteresis of (a) the energy E and (b) S(n,0), as 
V2 is increased and decreased showing the striped solid to 
superfluid transition to be first order. 



rying out the appropriate finite size scaling. What this 
analysis will show is that while a quite reasonable data 
collapse can be achieved, the requisite dynamic critical 
exponent is anomalously small. This will lead us to per- 
form a more careful analysis which will reveal the true, 
and more subtle, critical behavior. 

In Fig. 4 we show p s versus V2 for system sizes L = 
6, 8, 10, 12. In a quantum phase transition, the finite size 
scaling function not only depends on the appropriately 
scaled distance to the critical point, but also on the ratio 
of the lattice dimensions in space and imaginary time. It 
is standard to assume the following form for the super- 
fluid densitytl, 



1 / y 9 _ yc 
oc —F [ 2 : 2 , 0/ L z 



(3) 



One approach is to simulate several sets of lattice sizes, 
each set with a different space/imaginary time aspect 
ratio associated with different guesses for z. Instead, we 
first use another approach which is to choose the inverse 
temperature large enough ((3 = 20 usually suffices) so 
that the second argument in F is a constant as L is varied. 
With this choice, p s L z must then be independent of L 
at the critical point. The best intersection was found 
when z = 0.25, giving a critical interaction value V 2 C — 
0.765 (see Fig. 4). Replotting this figure using the scaled 
variable (V2 — V 2 C )F X I V we obtain very good data collapse 
shown in figure Fig. 5. This collapse yields the value of 
the correlation function exponent v = 0.36. 

Once these critical exponents were found, we redid the 
simulations but with the two parameter finite size scal- 
ing analysis mentioned above in which the temperature 
is varied with L, (3 oc L z , in order to keep the second 
argument in Eq. || constant. This reproduced the same 
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FIG. 4. p s as a function of V2 for L = 6, i 
intersection gives the critical value V 2 C = 0.765. 
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values found for z and v. 

Scaling analyses and data collapse of this level of qual- 
ity are rather commonly used to draw conclusions con- 
cerning the appropriate critical behavior and universality 
However, the value of the dynamical critical exponent, 
z = 0.25, is surprisingly small. All previous debateE3'E3 
had been whether z = 1 or z — 2. This leads us to re- 
examine the process and, in particular, check the valid- 
ity of the above finite size scaling analysis, by applying 
exactly the same assumptions and methodology to the 
transition at the Heisenberg point, where we know they 
should not hold. We obtain very similar scaling data, 
with the same values of z and v. 

The similarity between the behavior of the superfluid 
density near the Heisenberg point with that at finite V2 
suggests that the entire line of phase transitions sepa- 
rating the SF and (tt, 7r)-solid phases might be in the 
same class as the Heisenberg point. However, especially 
given the possibility of generating 'acceptable' finite size 
scaling plots despite the known critical behavior, that we 
have just demonstrated, this conjecture clearly requires 
careful testing, which we shall now describe. 

The Heisenberg point is very special in that the Hamil- 
tonian has an 0(3) symmetry which is explicitly broken 
everywhere else in the (V^Va) plane. In order for the 
transition line from the SF to the (tt, 7r)-solid to be in 
the same universality class, this symmetry must be dy- 
namically restored at the transition. To check this nu- 
merically, we first need to determine accurately a tran- 
sition point away from the special Heisenberg point. To 
this end, we fixed V\ = 3 and studied, for many val- 
ues of V2 and L, the behavior of S(tt,tt) and the con- 
densate, N(k = 0) = af(k = 0)a(k = 0) where a is the 
Fourier transform of the destruction operator. The re- 
sults for N(k = 0) are shown in Fig. 6 and indicate 
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FIG. 5. Same as Fig. 4 but using the scaling variables 
PsL 25 and \V 2 - V£\ y ". This yields the values z = 0.25 and 
v = 0.36. 



that the transition happens at 0.760 < V2 < 0.762 for 
Vi = 3. Furthermore, a substantial jump in N(k) is indi- 
cated over this narrow window of parameters. This very 
abrupt transition argues for a discontinuous transition 
like that at the Heisenberg point. 

A sensitive test of this suggestion is that if the 0(3) 
symmetry is indeed restored at the critical V2, the tran- 
sition can only take place at zero temperature. To verify 
this we did simulations at finite temperature to deter- 
mine the transition temperatures as a function of V2 ■ For 
V2 < 0.761, the transition is from a normal bosonic liquid 
to the (tt, tt) solid, and is expected to be in the 2d Ising 
universality class. The transition temperatures were de- 
termined from crossing points of the 4-th order Binder 
cumulant ratios 1 — (S(ir, tt) 2 ) /3(S(n, tt)) 2 for different 
system sizes L, and plotted in Figure 7. For V 2 > 0.761, 
the transition is from a normal bosonic liquid to a su- 
perfluid, and is expected to be in the 2d XY universal- 
ity class. The critical temperaturs Tkt of this phase 
were determined from the universal jump of the super- 
fluid density p s at the critical temperature. In Fig. 7 we 
see very clearly that the transition temperature plunges 
to zero as we approach the critical V2. In addition, the 
way the critical temperature drops to zero is well fitted 
by — l/ln(V2 — V2c), just as is the case at the Heisenberg 
point. 

It seems very likely, therefore, that in the ground state, 
T = 0, the transition at half filling from the checkerboard 
solid to the superfluid proceeds via a dynamical restora- 
tion of symmetry and is in the same class as the special 
Heisenberg point at Vi = 2, V2 = 0. However, our simula- 
tions were done at finite, albeit very low, temperatures. 
Consequently, for the moment, we cannot exclude the 
possibility that at very low temperatures the transition 
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FIG. 6. The number of bosons in the zero momentum 
mode, the condensate, as a function of 1/L for different next 
near neighbor, V2, repulsion values. f3 = 2L. 
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FIG. 7. The transition temperature, T, versus Vi- At the 
critical V%{tt 0.761) the transition temperature drops to zero. 
This suggests possible restoration of the 0(3) symmetry. 



is first order. These issues and the finite temperature 
phase transition are currently being studied further r3 

It is interesting to compare this with what happens 
in other systems. In the 2-d extended fermion Hubbard 
model, there is a similar competition between an anti- 
ferromagnetic ordered phase and a charge density wave 
phase. In the case of fermions, the entire zero tempera- 
ture phase diagram in the (£/, Vi) (onsite and near neigh- 
bor repulsion) plane at half-filling is believed to be in- 
sulating. There is no metallic or superconducting re- 
gion. This is thought to arise in part because of the 
peculiar nature of the two-dimensional Fermi surface, 
which has perfect nesting and a logarithmic divergence of 
the density of states at half-filling. Both these features 
act to enhance the tendency for diagonal long range or- 
der. However, little has been pinned down precisely as 
to the nature of the phase transitions in two dimensions. 
In one dimension, the question has been addressed, and 
the transition between the two insulating states has been 
argued to be second order at weak coupling and first 
order patr-Sitrong coupling, with an intervening tricritical 
point E3€3 Recent work has called this into question, and 
suggested that, an intervening bond-ordered-wave phase 
is important .EJ 

In the 1-d soft core boson-Hubbard model, the ground 
state phase diagram in the U (on-site repulsion) and |4i 
(near neighbor repulsion) was also studied cxtcnsivclyllZI 
at filling n = 1. As in the case of 2-d reported here, 
there is a weak-coupling superfluid phase which is sup- 
planted by ordered, insulating Mott and charge density 
wave phases at large U and V\ respectively. A study 
of hysteresis loops and the free energy barriers indicated 
that the transition out of the cdw phase to the superfluid 
is second order at weak coupling, and that the transi- 
tion from cdw to Mott phase at strong coupling is first 
order. We see that there are similarities with the 2-d 
hardcore case reported here, but there also differences, 
namely here there could be a dynamically restored 0(3) 
symmetry which was explicitly broken in the hamilto- 



V. DOPED SYSTEM 

Now we examine the phase transitions and various 
phases when the system is doped away from half fill- 
ing. The hard-core boson-Hubbard Hamiltonian, Eq. 1 
has particle-hole symmetry. That is, the transformation 
a>i —> a. maps n\ — > 1 — m, interchanging occupied and 
empty sites, but leaves the Hamiltonian unchanged apart 
from trivial constants. Therefore it is sufficient to do the 
simulations below half-filling and use this symmetry to 
calculate physical quantities above. 
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FIG. 8. The structure factor, S(tv, tt), and superfluid den- FIG. 9. The particle density, p, as a function of the calcu- 

sity, p s , versus the particle density, p. For 0.4 < p < 0.5, both lated chemical potential, p. The slope is the compressibility, 

p a and S(tt,tt) are finite. k = dp/ dp. 



A. Evolution of Checkerboard Solid 

To explore what happens to the checkerboard solid 
when the system is doped, we performed a series of sim- 
ulations at various values of V\ as the number of bosons 
in the system is lowered from half filling. We did this 
mostly with the canonical dual QMC algorithm where 
the number of bosons is fixed and the chemical poten- 
tial is calculated from the energy to add a particle to the 
system, 

fi = E(N h 

O.SOll + 1) - E(N h 

os on ) (4) 

In Fig. 8 we show a typical result for V\ = 3. At low 
densities (in this case when p < 0.4) we see that p s is 
finite while S(tt, tt) is small and decreasing as we move 
away from p = 0.5. In addition, finite size studies show 
that p s is essentially unchanging while S(tt, tt) — > for 
a fixed p < 0.4 as L grows. Therefore, this corresponds 
to a superfluid phase. The structure factor reaches its 
maximum at p = 0.5 while p s is zero there. This is the 
checkerboard solid discussed in the previous section. 

Between the superfluid phase and the half filling 
checkerboard solid, i.e. for 0.4 < p < 0.5, Fig. 8 shows 
that both the structure factor and superfluid density are 
non- vanishing. Furthermore, this is not a finite size ef- 
fect: For larger systems, p s maintains its value while 
S(tt, tt) diverges with L 2 as it should in the case of long 
range density wave order. This, therefore, is a candidate 
for a checkerboard supersolid phase. 

To verify this possibility, and check the thermody- 
namic stability of the supersolid phase we show in Fig. 9 
the density, p, as a function of the calculated chemi- 
cal potential, p. We see that for all the density values 
where Fig. 8 shows a supersolid, i.e. 0.4 < p < 0.5, the 



curve in Fig. 9 has negative slope and therefore negative 
compressibility, n = dp/dp. Consequently, the apparent 
checkerboard supersolid phase is not stable thermody- 
namically and undergoes phase separation into a mixture 
of checkerboard solid and superfluid. This same behav- 
ior had previously been established for the magnetization 
process of the spin-1/2 XXZ model on smaller lattices.til 

To establish this phase separation further, we simu- 
lated the system in the grand canonical ensemble where 
p, is the input parameter and p is calculated. If the sys- 
tem undergoes phase separation, as shown in Fig. 9, then, 
for the corresponding value of p, a histogram of the den- 
sity should show two peaks, one at p = 0.5 and the other 
at p < 0.5. This is indeed what happens as shown in 
Fig. 10 for an 8 x 8 system at V\ = 2.86. The simulation 
is done for several values of the chemical potential. The 
phase transition takes place for the p, value with equal 
peaks. We verified that the peak separation does not 
change when the system size is increased. 

By repeating the simulations that led to Fig. 9 for var- 
ious values of V\, we map out the phase diagram in the 
(t/Vi, n/Vi) plane for V2 = 0. The transitions between 
the superfluid and the (tt, 7r)-solid phases are first or- 
der except at half filling which is the special Heisenberg 
point. This is shown in Fig. 11. This phase diagram is 
in a gree ment with the mean field/spin wave analysis in 
Ref. ||. 

B. Evolution of Striped Solid 

Now we investigate the effect of doping on the striped 
solid phase present at half filling. 

The top part of Fig. 12 shows the structure factors, 
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FIG. 10. Histogram of the particle density as the chem- 
ical potential, /i is changed. The double peaks show phase 
separation. 
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FIG. 11. The phase diagram for V 2 = 0. The solid line 
shows the continuous transition to the Mott phase at full fill- 
ing, the dashed line shows the discontinuous first order tran- 
sitions from the superfluid to the checkerboard solid at half 
filling. The tip of the lobe, p — 0.5, is the Heisenberg point. 
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FIG. 12. Top: 5(0,71-) (circles) and S(ir, 0) (triangles). 
Bottom: (W%) (circles) and (Wy) (triangles). The larger 
(W 2 ) (circles) is parallel to the stripes, the lower is trans- 
verse. The system is 8 x 8, Vi = 0, V 2 = 5, j3 = 6 



S(n, 0) and 5(0, tt), for Vi = 0, V 2 = 5. For p < 0.3 we see 
that the system is isotropic, S(tt, 0) = S(0, 7r) and vanish- 
ing. For p > 0.3 the symmetry is broken and one of the 
two vanishes while the other is large (diverges with the 
system size). That signals the formation of stripes along 
the x or y directions. It is remarkable that the stripes 
start forming at such small densities. The figure also 
shows that at the same densities where the stripes form, 
the superfluid density is no longer isotropic, p% ^ p v s . 
The superfluid densities in the x and y directions are de- 
fined by: p% = (W*)/4tp and p y s = {W%)/At(3, where 
W x {W v ) is the winding number in the x (y) direction. 
In addition, the figure shows clearly that the supcrlfuid 
density along the stripes is larger than transverse to the 
stripes. Nonetheless, the superfluid density in the trans- 
verse direction does not vanish. Therefore, once again 
we apparently have a phase which is both superfluid and 
solid, thus another candidate for the supersolid phase. 

Again we check the thermodynamic stability of the su- 
persolid phase by calculating p as a function of p. This is 
shown in Fig. 13. We see that the compressibility (slope) 
never becomes negative indicating that phase separation 
is absent. In addition, we see that the compressibility in- 
creases sharply at p = 0.3 which is the density at which 
stripes form (see Fig. 12). This indicates that, contrary 
to the checkerboard case, the striped supersolid phase is 
indeed thermodynamically stable and has a higher com- 
pressibility than the superfluid phase. 

It is worth emphasizing that the striped supersolid is 
not merely a one dimensional superfluid phase along the 
channels created by the stripes. If this were the case, 
the superfluid density transverse to the stripes would be 
vanishingly small, which it is not. In addition, it was 
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FIG. 13. Particle density, p, versus chemical potential, 
fi. There is a sharp increase in the compressibility as p is 
increased when the system goes from the superfluid to the 
supersolid phase. 



shown in Ref. ^Ojthat (a^ (r)a(f ')} is finite as \r— f'\ — » oo 
transverse to the stripes in the supersolid phase. 

By repeating the simulations that led to Figs. 12 and 
13, we map out the phase diagram in the (t/V2, (J./V2) 
plane. This is shown in Fig. 14. The narrow regions 
sandwiched between SF and (71" , 0) solid phases are the 
stable supersolid phases. 

As can be seen in Fig. 12, the transition from the SF 
to the supersolid phase appears to be first order. The 
transition from the supersolid to the striped solid phase 
is continuous as is seen from the behavior of the super- 
fluid density as half filling is approached. We see from 
the lower part of Fig. 12 that both branches of the super- 
fluid density, parallel and transverse to the stripes, go to 
zero smoothly as half filling is approached. In fact, both 
branches behave like \p — 1/2 1, indicating a second order 
transition with a unit exponent. 



VI. CONCLUSIONS 
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t/v 2 

FIG. 14. The phase diagram for Vi = 0. The narrow 
regions sandwiched between SF and (it, 0) solid phases are 
the stable supersolid phases. 



that, for the bosonic Hubbard model at half filling, there 
is no supersolid phase, between the checkerboard solid 
and gnporflitiH phasealj, unlike what is observed in other 
modelsB3o. The details of the finite tepaperature phase 
diagram are currently being worked outo. 

In addition, by examining the density histograms from 
a grand canonical algorithm, we verified again that what 
was previously thought to be a checkerboard supersolid 
phase, is in fact a phase separated mixture of superfluid 
and solid regions. 

As for the striped phases, we showed, from the hystere- 
sis of the energy and structure factor, that the superfluid 
to striped solid transition, at half filling, is first order. 
The transition from the superfluid phase to the striped 
supersolid phase (away from half filling) also appears to 
be first order, in disagreement with Refs. 31 p2. We also 



showed, from the behavior of the superfluid density, that 
the transition from the striped supersolid phase to the 
striped solid phase is second order with the superfluid 
density vanishing as p s ~ \p — 1/2 1. 



The boson-Hubbard model exhibits many fascinating 
quantum phases and phase transitions. In this paper, we 
have shown that the detailed critical behavior at those 
transitions, that is, both the order of the transitions and 
the critical exponents, can now be determined with re- 
cently developed Quantum Monte Carlo algorithms. Our 
principal conclusion is that, although a first order tran- 
sition cannot be categorically ruled out at the moment, 
the superfluid to checkerboard solid transition at half- 
filling appears to proceed via a dynamical restoration of 
the explicitly broken 0(3) symmetry, and is therefore in 
the same class as the Heisenberg point. This confirms 



Acknowledgments 
We acknowledge useful conversations with E.W. Fire. 
This work was supported by NSF-DMR-9985978. We 
also thank ETH-Zurich and HLRS (Stuttgart) for very 
generous grants of computer time. 







1 M.P.A. Fisher, P.B. Weichman, G. Grinstein, and D.S. 
Fisher, Phys. Rev. B40, 546 (1989). 

2 M. Cha, M.P.A. Fisher, S.M. Girvin, M. Wallin, and A.P 
Young, Phys. Rev. B44, 6883 (1991). 

3 E.S. Sorensen, M. Wallin, S.M. Girvin, and A.P. Young, 
Phys. Rev. Lett. 69, 828 (1992). 

4 G.G. Batrouni, B. Larson, R.T. Scalettar, J. Tobochnik, 
and J. Wang, Phys. Rev. B48, 9628 (1993). 

5 K.J. Runge, Phys. Rev. B45, 13136 (1992). 

6 V.F. Elesin, V.A. Kashurnikov, and L.A. Openov, JETP 
Lett. 60, 177 (1994). 

7 A.P. Kampf and G.T. Zimanyi, Phys. Rev. B47, 279 
(1993). 

8 R.T. Scalettar, G.G. Batrouni, and G.T. Zimanyi, 
Phys. Rev. Lett. 66, 3144 (1991). 

9 G.T. Zimanyi, P.A. Crowell, R.T. Scalettar, and G.G. Ba- 
trouni, Phys. Rev. B50, 6515 (1994). 

10 K.G. Singh and D.S. Rokhsar, Phys. Rev. B46, 3002 
(1992). 

11 M. Kohno and M. Takahashi, Phys. Rev. B56, 3212 (1997). 

12 G.G. Batrouni, R.T. Scalettar, G.T. Zimanyi and A.P. 
Kampf, Phys. Rev. Lett. 74 (1995) 2527; R.T. Scalettar, 
G.G. Batrouni, A.P. Kampf and G.T. Zimanyi, Phys. Rev. 
B51 (1995) 8467. 

13 G.G. Batrouni and R.T. Scalettar, Phys. Rev. Lett. 84 
1599 (2000). 

14 T.D. Kuhner, S.R. White, and H. Monien, Phys. Rev. B61, 
12474 (2000). 

15 M.C. Cha and S.M. Girvin, Phys. Rev. B49, 9794 (1994). 

16 G.G. Batrouni, R.T. Scalettar, and G.T. Zimanyi, 
Phys. Rev. Lett. 65, 176 (1990). 

17 P. Niyaz, R.T. Scalettar, CY. Fong, and G.G. Ba- 
trouni, Phys. Rev. B44, 7143 (1991); G.G. Batrouni 
and R.T. Scalettar, Phys. Rev. B46, 9051 (1992); and 
P. Niyaz, R.T. Scalettar, CY. Fong, and G.G. Batrouni, 
Phys. Rev. B50, 362 (1994). 

18 W. Krauth and N. Trivedi, Europhys. Lett. 14, 627 (1991). 

19 G.G. Batrouni and H. Mabilat, Comp. Phys. Comm. 121- 
122 468 (1999). 

20 F. Hebert, G.G. Batrouni, and H. Mabilat, Phys. Rev. B61 
10725 (2000). 

21 G.G. Batrouni, Nucl. Phys. B208 467 (1982); G.G. Ba- 
trouni and M. B. Halpern, Phys. Rev. D30 1775 (1984). 

22 A. Sandvik, Phys. Rev. B59, R14157 (1999). 



2:s 



A. Dorneich and M. Troyer, unpublished. 



24 H.G. Evertz, G. Lana, and M. Marcu, Phys. Rev. Lett. 70, 
875 (1993); U.-J. Wiese and H.-P. Ying, Phys. Lett. A 168, 
143 (1992); Z. Phys. B 93 (1994) 147; N. Kawashima, J.E. 
Gubernatis, and H.G. Evertz, Phys. Rev. B 50 136 (1994); 
N. Kawashima and J.E. Gubernatis, Phys. Rev. Lett. 73, 
1295 (1994), B. B. Beard and U.-J. Wiese, Phys. Rev. Lett. 
77, 5130 (1996). 

25 Some additional studies of the boson-Hubbard model with 
disorder include: N. Hatano, Int. J. Mod. Phys. C7, 449 
(1996); N. Hatano, J. Phys. Soc. Japan 64, 1529 (1995); 
J. Kisker and H. Rieger, Phys. Rev. B55, 11981 (1997); J. 
Kisker and H. Rieger, Physica A246, 348 (1997); 

26 E.L. Pollock and D.M. Ceperley, Phys. Rev. B30, 2455 
(1984) and Phys. Rev. B36, 8343 (1987). D.M. Ceperley 
and E.L. Pollock, Phys. Rev. Lett. 56, 351 (1986). 



27 A. van Otterlo, K.-H. Wagenblast, Phys. Rev. Lett. 72, 
3598 (1994), A. van Otterlo et al, Phys. Rev. B52, 16176 
(1995) 

28 J.E. Hirsch, Phys. Rev. Lett. 53 2327 (1984). Phys. Rev. 
B31 (1991) 6022. 

29 J.W. Cannon, R.T. Scalettar, and E. Fradkin, Phys. Rev. 
B44 (1991) 5995. 

30 P. Sengupta, A. Sandvik, and D.K. Campbell, cond- 
mat//0102141. 

31 E. Frey and L. Balents, Phys. Rev. B55, 1050 (1997). 

32 C. Pich and E. Frey, Phys. Rev. B57, 13712 (1998). 

33 G. Schmid, M. Troyer and A. Dorneich, unpublished. 

34 E. Roddick and D. Stroud, Phys. Rev. B51, 8672 (1995). 



10 



